#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkMath.h>
#include <vtkPolyDataReader.h>
#include <vtkSmoothPolyDataFilter.h>
#include "vtkScalarBarActor.h"
#include <vtkLookupTable.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkLookupTable.h>

#include <vtkCenterOfMass.h>
#include <vtkAssembly.h>
#include <vtkSphereSource.h>
#include "vtkCellArray.h"
#include "vtkImplicitModeller.h"
#include "vtkContourFilter.h"
#include "vtkTransform.h"
#include "vtkTransformFilter.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkWarpTo.h"
#include "vtkDataSetMapper.h"
#include "vtkRotationalExtrusionFilter.h"
#include "vtkCamera.h"
#include <vtkOrientationMarkerWidget.h>
#include "vtkAnnotatedCubeActor.h"
#include "vtkAxesActor.h"
#include "vtkTransform.h"
#include "vtkTextProperty.h"
#include "vtkCaptionActor2D.h"
#include "vtkPropAssembly.h"
#include "vtkLODActor.h"


#define VTK_CREATE(type, var) \
  vtkSmartPointer<type> var = vtkSmartPointer<type>::New()

#define USAGE "Usage: CurvatureMapVisualization 1 2 3 4\n"\
			  "where\n"\
			  "\t1. CurvatureMap1 - Generated CurvatureMap description of the object1 (Diastole) \n"\
			  "\t2. CurvatureMap2 - Generated CurvatureMap description of the object2 (Systole)\n"\
			  "\t3. width - specify the window width in pixels \n"\
			  "\t4. height - specify the window height in pixels'\n"\
              "\t5. inSliceSpacing - Optional: inSliceSpacing = X-spacing = Y-spacing'\n"\
              "\t6. x1 position - Optional: specify the x pixel-position of reference point on object1'\n"\
              "\t7. y1 position - Optional: specify the y pixel-position of reference point on object1'\n"\
              "\t8. x2 position - Optional: specify the x pixel-position of reference point on object2'\n"\
              "\t9. y2 position - Optional: specify the y pixel-position of reference point on object2'\n"

// ./VisualizeCurvatureMap ./curvature9.vtk  ./curvature18.vtk 1024 512 1.171875 65 125 65 125;
// ./VisualizeCurvatureMap ./SPHARM1.vtk  ./SPHARM10.vtk 1024 512
// ./VisualizeCurvatureMap ./SPHARM1.vtk  ./SPHARM10.vtk 1024 512 1.171875 65 125 65 125;

void AttributeAdjustment(vtkPolyData* inputPolyData)
{

	int nTuples = vtkDoubleArray::SafeDownCast(	inputPolyData->GetPointData()->GetScalars() )->GetNumberOfTuples();
	std::cout << "nTuples: " << nTuples << std::endl;
	double* scalars = new double[nTuples];
	double mScalars = 0;
	for (int i = 0; i < nTuples; i++) // Calculate:  mean(x(:))
	{
		scalars[i] = vtkDoubleArray::SafeDownCast( inputPolyData->GetPointData()->GetScalars() )->GetValue(i);
		mScalars += scalars[i];
//std::cout << "scalars[i]: " << scalars[i] << std::endl;
	}
	mScalars /= nTuples;
	std::cout << "1--mScalars: " << mScalars << std::endl;

	for (int i = 0; i < nTuples; i++) // Update:  x = x-mean(x(:));
	{
		scalars[i] = scalars[i] - mScalars;
	}

	mScalars = 0;
	for (int i = 0; i < nTuples; i++) //ReCalculate:  mean(x(:))
	{
		mScalars += scalars[i];
	}
	mScalars /= nTuples;
	std::cout << "2--mScalars: " << mScalars << std::endl;

	double madScalars = 0;
	for (int i = 0; i < nTuples; i++) // Calculate:  mad(X(:)) = MEAN(ABS(X-MEAN(X))
	{
		madScalars += fabs(scalars[i] - mScalars);
	}
	madScalars /= nTuples;
	std::cout << "madScalars: " << madScalars << std::endl;

	for (int i = 0; i < nTuples; i++) // Update:  x = x/(2*mad(x(:)));
	{
//		scalars[i] = scalars[i] / (2 * madScalars);
//		if (scalars[i] > 2)
//			scalars[i] = 2;
//		if (scalars[i] < -2)
//			scalars[i] = -2;
		scalars[i] = scalars[i] / (2 * madScalars);
		if (scalars[i] > 1)
			scalars[i] = 1;
		if (scalars[i] < -1)
			scalars[i] = -1;
//		if (scalars[i] < 0)
//			scalars[i] = -scalars[i];
	}

	for (int i = 0; i < nTuples; i++) // Calculate:  mean(x(:))
	{
		vtkDoubleArray::SafeDownCast( inputPolyData->GetPointData()->GetScalars() )->SetValue(i,scalars[i]);
	}

}

int main(int argc, char *argv[])
{


	float inSliceSpacing=0;
	int x1Position=0, y1Position=0, x2Position=0, y2Position=0;

	if (argc > 9)
	{
	inSliceSpacing = atof(argv[5]);
	x1Position = atoi(argv[6]);
	y1Position = atoi(argv[7]);
	x2Position = atoi(argv[8]);
	y2Position = atoi(argv[9]);
	}

	vtkPolyDataReader *reader1 = vtkPolyDataReader::New();
	reader1->SetFileName(argv[1]);
	reader1->Update();

	vtkPolyDataReader *reader2 = vtkPolyDataReader::New();
	reader2->SetFileName(argv[2]);
	reader2->Update();

	int width, height;
	width = atoi(argv[3]);
	height = atoi(argv[4]);

	vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
//	lookupTable->SetTableRange(-1.0, 1.0);
	lookupTable->SetTableRange(0, 0.15);
	lookupTable->SetNumberOfColors(512);
	lookupTable->SetNumberOfTableValues(512);
	lookupTable->Build();
    for (int i = 0; i < 512; i++)
    {
    	double 	rgba[4];
    	lookupTable->GetTableValue(i, rgba);
        lookupTable->SetTableValue(i, rgba[2], rgba[1], rgba[0], rgba[3]);
    }



	 // Object 1
	vtkSmartPointer < vtkPolyData > inputPolyData1 = reader1->GetOutput();
//	AttributeAdjustment(inputPolyData1);

	// Top slice location
	double bounds[6];
	inputPolyData1->GetBounds(bounds);
	double polyData1Zmax = bounds[5];
//	  std::cout  << "xmin: " << bounds[0] << " "
//	             << "xmax: " << bounds[1] << std::endl
//	             << "ymin: " << bounds[2] << " "
//	             << "ymax: " << bounds[3] << std::endl
//	             << "zmin: " << bounds[4] << " "
//	             << "zmax: " << bounds[5] << std::endl;

	// Compute the center of mass
	vtkSmartPointer < vtkCenterOfMass > centerOfMassFilter = vtkSmartPointer< vtkCenterOfMass > ::New();

#if VTK_MAJOR_VERSION <= 5
	centerOfMassFilter->SetInput(inputPolyData1);
#else
	centerOfMassFilter->SetInputData(inputPolyData1);
#endif
	centerOfMassFilter->SetUseScalarsAsWeights(false);
	centerOfMassFilter->Update();

	double center1[3];
	centerOfMassFilter->GetCenter(center1);
	std::cout << "Mass center of Object1 is: " << center1[0] << " " << center1[1]<< " " << center1[2] << std::endl;


	// Create input mapper
	vtkSmartPointer < vtkPolyDataMapper > inputMapper1 = vtkSmartPointer< vtkPolyDataMapper > ::New();
#if VTK_MAJOR_VERSION <= 5
  inputMapper1->SetInput(inputPolyData1);
#else
  inputMapper1->SetInputData(inputPolyData1);
#endif

//    double scalarRange[2];
//    inputMapper1->GetScalarRange(scalarRange);
//    std::cout << "scalarRange is: " << scalarRange[0] << ", " << scalarRange[1]<< "; "<< std::endl;

    inputMapper1->SetScalarRange(0,0.15);
//    inputMapper1->SetScalarRange(-1,1);
	inputMapper1->SetLookupTable(lookupTable);

//	inputMapper1->GetScalarRange(scalarRange);
//	std::cout << "scalarRange is: " << scalarRange[0] << ", " << scalarRange[1]<< "; "<< std::endl;

	// Create input actor
	vtkSmartPointer < vtkActor > inputActor1 = vtkSmartPointer < vtkActor > ::New();
	inputActor1->SetMapper(inputMapper1);
	inputActor1->SetPosition(-center1[0],-center1[1],-center1[2]);


	// Object 2
	vtkSmartPointer < vtkPolyData > inputPolyData2 = reader2->GetOutput();
//	AttributeAdjustment(inputPolyData2);

	// Top slice location
	inputPolyData2->GetBounds(bounds);
	double polyData2Zmax = bounds[5];


	// Compute the center of mass
#if VTK_MAJOR_VERSION <= 5
	centerOfMassFilter->SetInput(inputPolyData2);
#else
	centerOfMassFilter->SetInputData(inputPolyData2);
#endif
	centerOfMassFilter->SetUseScalarsAsWeights(false);
	centerOfMassFilter->Update();

	double center2[3];
	centerOfMassFilter->GetCenter(center2);
	std::cout << "Mass center of Object2 is: " << center2[0] << " " << center2[1]<< " " << center2[2] << std::endl;


	// Create input mapper
	vtkSmartPointer < vtkPolyDataMapper > inputMapper2 = vtkSmartPointer < vtkPolyDataMapper > ::New();
#if VTK_MAJOR_VERSION <= 5
  inputMapper2->SetInput(inputPolyData2);
#else
  inputMapper2->SetInputData(inputPolyData2);
#endif

//    inputMapper2->SetScalarRange(-1,1);
    inputMapper2->SetScalarRange(0, 0.15);
	inputMapper2->SetLookupTable(lookupTable);

	// Create input actor
	vtkSmartPointer < vtkActor > inputActor2 = vtkSmartPointer < vtkActor > ::New();
	inputActor2->SetMapper(inputMapper2);
	inputActor2->SetPosition(-center2[0], -center2[1], -center2[2]);



//	  inputActor2->SetOrigin(center[0],center[1],center[2]);
	/*
	double rrange[2];
	inputPolyData->GetPointData()->GetScalars()->GetRange(rrange);
	std::cout << "rrange: " << rrange[0] << " , " << rrange[1] << std::endl;

	double range[2];
	inputMapper->GetScalarRange(range);
	std::cout << "range: " << range[0] << " , " << range[1] << std::endl;

    //inputMapper->SetScalarRange(rrange[0], 0);
	*/


	// Coordinate system
	vtkPolyData *pd2 = vtkPolyData::New();
	vtkCellArray *ca2 = vtkCellArray::New();
	vtkPoints *pts2 = vtkPoints::New();
	pts2->InsertNextPoint(0, 1, 0);
	pts2->InsertNextPoint(8, 1, 0);
	pts2->InsertNextPoint(8, 2, 0);
	pts2->InsertNextPoint(10, 0.01, 0);
	ca2->InsertNextCell(4);
	ca2->InsertCellPoint(0);
	ca2->InsertCellPoint(1);
	ca2->InsertCellPoint(2);
	ca2->InsertCellPoint(3);
	pd2->SetPoints(pts2);
	pd2->SetLines(ca2);

	vtkTransform *arrowT2 = vtkTransform::New();
	arrowT2->Scale(1, 0.6, 1);
	arrowT2->RotateY(90);

	vtkTransformPolyDataFilter *arrowTF2 = vtkTransformPolyDataFilter::New();
	arrowTF2->SetInput(pd2);
	arrowTF2->SetTransform(arrowT2);

	vtkRotationalExtrusionFilter *arrowREF = vtkRotationalExtrusionFilter::New();
	arrowREF->SetInput(arrowTF2->GetOutput());
	arrowREF->CappingOff();
	arrowREF->SetResolution(30);

	vtkPolyDataMapper *spikeMapper = vtkPolyDataMapper::New();
	spikeMapper->SetInput(arrowREF->GetOutput());

	vtkLODActor *a1Actor = vtkLODActor::New();
	a1Actor->SetMapper(spikeMapper);
	a1Actor->SetScale(1, 1, 5);
	a1Actor->RotateY(180);
	a1Actor->SetPosition(0, 0, 0);
	a1Actor->GetProperty()->SetColor(1,1, 1);
	a1Actor->GetProperty()->SetSpecularColor(1,1, 1);
	a1Actor->GetProperty()->SetSpecular(0.3);
	a1Actor->GetProperty()->SetAmbient(0.2);
	a1Actor->GetProperty()->SetDiffuse(0.8);
	a1Actor->GetProperty()->SetSpecularPower(20);

	vtkLODActor *a2Actor = vtkLODActor::New();
	a2Actor->SetMapper(spikeMapper);
	a2Actor->SetScale(1, 1, 4);
	a2Actor->RotateY(-90);
	a2Actor->SetPosition(0, 0, 0);
	a2Actor->GetProperty()->SetColor(1,1, 1);
	a2Actor->GetProperty()->SetSpecularColor(1,1, 1);
	a2Actor->GetProperty()->SetSpecular(0.3);
	a2Actor->GetProperty()->SetAmbient(0.2);
	a2Actor->GetProperty()->SetDiffuse(0.8);
	a2Actor->GetProperty()->SetSpecularPower(20);

	vtkLODActor *a3Actor = vtkLODActor::New();
	a3Actor->SetMapper(spikeMapper);
	a3Actor->SetScale(1, 1, 4);
	a3Actor->RotateZ(90);
	a3Actor->RotateY(-90);
	a3Actor->SetPosition(0, 1, 0);
	a3Actor->GetProperty()->SetColor(1,1, 1);
	a3Actor->GetProperty()->SetSpecularColor(1,1, 1);
	a3Actor->GetProperty()->SetSpecular(0.3);
	a3Actor->GetProperty()->SetAmbient(0.2);
	a3Actor->GetProperty()->SetDiffuse(0.8);
	a3Actor->GetProperty()->SetSpecularPower(20);

	// reference point
	vtkSphereSource *ss = vtkSphereSource::New();
	ss->SetRadius(2);
	vtkPolyDataMapper *rpMapper = vtkPolyDataMapper::New();
	rpMapper->SetInput(ss->GetOutput());

	vtkLODActor *rpActor1 = vtkLODActor::New();
	rpActor1->SetMapper(rpMapper);
	rpActor1->SetPosition(-x1Position*inSliceSpacing-center1[0], -y1Position*inSliceSpacing-center1[1], 5+polyData1Zmax-center1[2]);
	rpActor1->GetProperty()->SetSpecularColor(1,1,0);
	rpActor1->GetProperty()->SetSpecular(0.3);
	rpActor1->GetProperty()->SetAmbient(0.2);
	rpActor1->GetProperty()->SetDiffuse(0.8);
	rpActor1->GetProperty()->SetSpecularPower(20);

	vtkLODActor *rpActor2 = vtkLODActor::New();
	rpActor2->SetMapper(rpMapper);
	rpActor2->SetPosition(-x2Position*inSliceSpacing-center2[0], -y2Position*inSliceSpacing-center2[1], 5+polyData2Zmax-center2[2]);
	rpActor2->GetProperty()->SetSpecularColor(1,1,0);
	rpActor2->GetProperty()->SetSpecular(0.3);
	rpActor2->GetProperty()->SetAmbient(0.2);
	rpActor2->GetProperty()->SetDiffuse(0.8);
	rpActor2->GetProperty()->SetSpecularPower(20);


	// combine orientation markers into one with an assembly
	vtkSmartPointer < vtkPropAssembly > assembly1 = vtkSmartPointer	< vtkPropAssembly > ::New();
	assembly1->AddPart(a1Actor);
	assembly1->AddPart(a2Actor);
	assembly1->AddPart(a3Actor);
	if (argc > 9) assembly1->AddPart(rpActor1);
	assembly1->AddPart(inputActor1);

	vtkSmartPointer < vtkPropAssembly > assembly2 = vtkSmartPointer	< vtkPropAssembly > ::New();
	assembly2->AddPart(a1Actor);
	assembly2->AddPart(a2Actor);
	assembly2->AddPart(a3Actor);
	if (argc > 9) assembly2->AddPart(rpActor2);
	assembly2->AddPart(inputActor2);




	VTK_CREATE(vtkScalarBarActor, scalarBar);
	scalarBar->SetLookupTable(lookupTable);
	scalarBar->SetTitle("Curvature");
	scalarBar->GetPositionCoordinate()->SetCoordinateSystemToNormalizedViewport();
	scalarBar->GetPositionCoordinate()->SetValue(0.1, 0.05);
	scalarBar->SetOrientationToVertical();
	scalarBar->SetWidth(0.8);
	scalarBar->SetHeight(0.9);
	scalarBar->SetLabelFormat("%-#6.3f");

	// Create renderers
	vtkSmartPointer < vtkRenderer > renderer1 = vtkSmartPointer < vtkRenderer> ::New();
	renderer1->AddActor(assembly1);
	renderer1->SetViewport(0.45, 0, 0.9, 1.0);
	renderer1->SetBackground(.1, .2, .3);

	VTK_CREATE(vtkRenderer, renderer2);
	renderer2->AddActor(assembly2);
	renderer2->SetViewport(0, 0, 0.45, 1.0);
	renderer2->SetBackground(.1, .2, .3);

	VTK_CREATE(vtkRenderer, renderer3);
	renderer3->AddActor2D(scalarBar);
	renderer3->SetViewport(0.9, 0, 1.0, 1.0);
	renderer3->SetBackground(0.3, 0.3, 0.3);



	//Add renderer to renderwindow and render
	vtkSmartPointer < vtkRenderWindow > renderWindow = vtkSmartPointer< vtkRenderWindow > ::New();
	renderWindow->AddRenderer(renderer1);
	renderWindow->AddRenderer(renderer2);
	renderWindow->AddRenderer(renderer3);
	renderWindow->SetSize(width, height);

	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
	renderWindowInteractor->SetRenderWindow(renderWindow);


    // Reposition Camera
	vtkCamera* cam1 = (renderer1->GetActiveCamera());
//	cam1->Azimuth(50);
	cam1->Elevation(-60);
	renderer1->ResetCamera();

	vtkCamera* cam2 = (renderer2->GetActiveCamera());
//	cam1->Azimuth(50);
	cam2->Elevation(-60);
	renderer2->ResetCamera();

	//Render and interact
	renderWindow->Render();
	renderWindowInteractor->Start();


	return EXIT_SUCCESS;
}
